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Abstract. We develop a model anisotropy best- 
fitting to the two-dimensional sky-map of multi-TeV 
galactic cosmic ray (GCR) intensity observed with 
the Tibet III air shower (AS) array. By incorpo- 
rating a pair of intensity excesses in the hydrogen 



deflection plane (HDP) suggested by Gurnett et al., 
together with the uni-directional and bi-directional 
flows for reproducing the observed global feature, 
this model successfully reproduces the observed sky- 
map including the "skewed" feature of the excess 
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intensity from the heliotail direction, whose physical 
origin has long remained unknown. These additional 
excesses are modeled by a pair of the northern and 
southern Gaussian distributions, each placed ~50° 
away from the heliotail direction. The amplitude of 
the southern excess is as large as ~0.2 %, more 
than twice the amplitude of the northern excess. 
This implies that the Tibet AS experiment discovered 
for the first time a clear evidence of the significant 
modulation of GCR intensity in the heliotail and the 
asymmetric heliosphere. 

Keywords: Best-fit model for the sidereal 
anisotropy, GCR modulation in the heliotail. 
Asymmetries of the heliosphere 

I. Introduction 

The Tibet III Air Shower Array experiment has been 
conducted at Yangbajing (90.522°E, 30.102°N; 4300 m 
above sea level) in Tibet, China. The array composed 
of 533 scintillation counters of 0.5 m^ each covers a 
detection area of 22,050 m^ achieving a trigger rate of 
~680 Hz. GCR events are selected for analyses, if all the 
following criteria are met: (i) any four-fold coincidence 
occurs in the counters with each recording more than 
0.8 particles in charge, (ii) the air shower core position 
is located in the array, (iii) the zenith angle of arrival 
direction is < 45°. With all these criteria, the array 
has the modal GCR energy of ~5 TeV. The angular 
resolution of the arrival direction of each shower is 
estimated to be ~0.9° from Monte Carlo simulations, 
and this was also verified by measuring the Moon's 
shadow in GCRs [IJ. In this paper, we analyze a total 
of 37 billion air shower events recorded in 1318.9 live 
days from November 1999 to October 2005. 

Fig. 1(a) shows the observed GCR intensity in 5° x 5° 
pixels in a color-coded format as a function of the right 
ascension (a) on the horizontal axis and the declination 
(6) on the vertical axis. For producing this figure, we first 
obtain the GCR intensity I^^^ in n-th right ascension 
and m-th declination pixel. We then normalize the 
average of I^^^ in each declination belt to unity and 
get the normalized model intensity i^f^ plotted in this 
figure. This sky-map covers 360° of a but covers only 
90° of 6 ranging from -15° to +75° due to the event 
selection criterion limiting zenith angles to < 45°. The 
map clearly shows a significant anisotropy, consisting of 
a --0.2 % excess at a 60° and J -10° and a --0.2 
% deficit at a ^ 180° and ^ ~ 0°, each observed with a 
statistical significance more than ten times the statistical 
error. There is also the "skewed" feature seen in the 
region of excess intensity as pointed by Q. In the next 
section, we develop a model anisotropy reproducing this 
observed global feature as well as the "skewed" feature. 

II. Analysis and result 

We develop a model anisotropy /n,m consisting of 
two components as 



where and respectively denote the global 

anisotropy (GA) and the additional excess (AE) intensity 
as described below. We first normalize the average of 
In,m in each declination belt to unity and get the 
normalized model intensity in,m = ^n,m + ^n,m' the 
same way as we did for the observed data to produce 
Fig. 1(a). By comparing in,m with z^^^ in Fig. 1(a), 
we obtain best-fit parameters minimizing the residual S 
defined, as 

Q — f^obs _ • \2 / _2 

where On.m is the statistical error of the intensity in 
(n, m) pixel. From this best-fitting, we excluded 16 
pixels containing the known and possible gamma ray 
sources. These pixels also include the "region A" re- 
ported by Milagro experiment |3|. 

The observed angular separation between the excess 
and the deficit in Fig. 1(a) is only ~ 120°, which is 
much smaller than 1 80° expected from a uni-directional 
flow (UDF) but significantly larger than 90° expected 
from a bi-directional flow (BDF). Only a combination of 
the uni-directional and bi-directional flows can achieve 
a reasonable fit to the global feature of the observed 
anisotropy. From this point of view, we perform a best- 
fit calculation to the observed global anisotropy (GA) 
with a model intensity expressed, as 

= ai±cosxi(n,m : 

+ai|| COSX2(^, ^ OL2,^2) 

+a2|| cos^ X2(n, m : 0^2, ^2) (3) 

where a\\_ and a^w are amplitudes of UDFs perpendic- 
ular and parallel to the BDF, respectively, a2\\ is the 
amplitude of the BDF, (ai, ^1) and (0^2, ^2) are respec- 
tively right ascensions and declinations of the reference 
axes of the perpendicular UDF and BDF and xi (X2) is 
the angle of the center of (n, m) pixel measured from the 
reference axis of the perpendicular UDF (BDF) |2|. Fig. 
1(b) displays derived from best-fitting to Fig. 1(a). 
While the model anisotropy in eq. (3) well reproduces 
the observed global feature in Fig. 1(a), it is obviously 
too simple for reproducing the observed feature of the 
excess intensity at around ol ~ 60° and J ~ —10°. 
This is demonstrated clearly in Fig. 1(c) displaying the 
significance of the residual anisotropy remaining after 
the subtraction of from Fig. 1(a). 

It is clear that the observed anisotropy contains an 
additional excess intensity along a plane which is almost 
perpendicular to the galactic plane. This additional ex- 
cess intensity was first reported as the "skewed" feature 
from the long-term observations of sub-TeV GCR inten- 
sity with underground muon detectors, but its physical 
origin has remained unknown |5||6|. It is also clear 
in Fig. 1(c) that the additional excess intensity extends 
along Gurnett's HDP displayed by a black curve |7|. We 
model this additional excess (AE) intensity by a pair of 
Gaussians placed in Gurnett's HDP, each centered away 
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Fig. 1. 2D-sky maps of the observed and reproduced GCR intensity. The panels display the normalized GCR intensity or significance in 
5° X 5° pixels in a color-coded format as a function of the right ascension on the horizontal axis and the declination on the vertical axis. In 
this figure, the average intensity in each declination belt is normalized to unity. These sky-maps cover 360° of the right ascension but cover 
only 90° of the declination due to the event selection criterion limiting zenith angles to < 45° (see text). The data in 16 pixels containing 
the known and possible gamma ray sources are excluded from the best-fit calculation and indicated by black color. The white curve indicates 
the galactic plane, while the black curve displays the HDP plane suggested by Gurnett et al. |7 |, which is calculated as a plane normal to the 
orientation of o; = 332.1° and 6 = 35.5°. In each panel, the heliotail direction (a = 75.9° and 6 = 17.4°) is indicated by a white solid circle 
in the HDP plane. Each panel displays, (a): the observed intensity (i^^^), (b): the best-fit component anisotropy (i^^,m) reproducing the global 
anisotropy (GA), (c): the significance of the residual anisotropy remaining after the subtraction of from i^^^ ( (i^^m ~ ^n,m)/<^ri,m), 
(d): the best-fit component anisotropy (in^) reproducing the additional excess (AE) intensity, (e): the total model anisotropy in,m best-fitting 
to (a), (f): the residual significance — in,m) /crn,m)- In panel (b), the open triangle with an attached character "F" indicates the LISMF 

orientation by Frisch (a = 300.9° and 6 = 32.2°) |12|, while the open diamond with "B" indicates the orientation of the best-fit BDF (see 
text and Table I). 

TABLE I 

Best- FIT parameters in 7^:^ and I^^^ in eqs. (3) and (4). Note that 62 in the upper table is derived according to our 

DEFINITION THAT THE REFERENCE AXES (ai,6i) AND (0:2,(^2) ARE PERPENDICULAR TO EACH OTHER [2|. 



aM%) ai||(%) a2||(%) 




^2(°) 


bl{%) b2{%) 


^ii(°) ^±n 


0.141 0.006 0.140 


37.5 37.5 


102.5 -28.9 


0.234 0.100 


25.0 10.0 52.5 



from the heliotail direction by an angle <l>, as 

= 1^1 exp( —2 ) 



+^2exp( —2 )\^M-^) (4) 

(f) 9 

where bi and 62 are amplitudes, and are widths 
parallel and perpendicular to the HDP respectively, 
0n,m is the "longitude" of the center of (n, m) pixel 
measured from the heliotail along the HDP and On,m 
is the "latitude" measured from the HDP. Fig. 1(d) 
displays the best-fit i^m^ while Fig. 1(e) displays the 
combined model anisotropy in,m best-fitting to i^f^ 
in Fig. 1(a). As seen in Fig. 1(f) showing the signif- 
icance of the residual anisotropy remaining after the 
subtraction of in,m from Fig. 1(a), this new model 
in,m successfully reproduces the observed anisotropy 
including the "skewed" feature of the excess intensity 
as well as the global feature. Eleven best-fit param- 
eters ( ai^,ai||,a2||,ai, ^1,^2, 61, 62, (70,(761,$) mini- 
mizing the residual 5* in eq. (2) are listed in Table 



I. The minimum S divided by the degree of freedom 
( 1 269=72 X 18-16-11) is 1.791. Note that five amplitudes 
(aix, tt2||7 ^2) are uniquely determined by the 
least- square technique for each set of remaining six 
angle parameters (ai, ^1 , 0^2 , cr^ , (751 , <l>). 

HI. Summary and discussion 

The Larmor radius of 5 TeV GCR protons in a 3 
/iG interstellar magnetic field is ~0.002 pc or 400 
AU, which is comparable to the scale of the heliosphere 
in the nose direction toward the upstream of the inter- 
stellar wind. This is one reason why the heliospheric 
modulation of the GCR intensity has been considered 
to be negligible in this energy region. The heliosphere, 
however, is also known to have a long heliotail extending 
over thousands of AU, much longer than of 5 TeV 
GCRs |8|. The GCR modulation in the heUotail still 
remains possible, although it is not fully understood yet. 

Recent observations also have suggested asymmetries 
of the heliosphere. Based on the deflection of the in- 
terstellar neutral hydrogen flow vector from the helium 
flow vector observed at 1 AU, Lallement et al. deduced 
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the Hydrogen Deflection Plane (HDP) containing the 
intersteUar wind velocity and the Local Interstellar Mag- 
netic Field (LISMF) and suggested a possible north- 
south asymmetry of the heliosphere due to the magnetic 
pressure of the LISMF |9|. From the observation of 2-3 
kHz radio emission from the outer heliosphere, Gurnett 
et al. |7| also deduced the HDP almost perpendicular 
to the galactic plane in a reasonable agreement with 
Lallement et al. Opher et al. analyzed the difference in 
the heliocentric distances to the solar wind termination 
shock observed by Voyager si and 2 at different he- 
liographic latitudes and longitudes. By comparing the 
difference with the Magneto-Hydro-Dymanic (MHD) 
simulations, they derived the north- south and east- west 
asymmetries of the heliosphere |10||11|. 

To achieve a good fit to the observed sky-map by 
Tibet AS experiment, the model anisotropy is required 
to include an additional excess intensity in the Gurnett 's 
HDP plane. This additional anisotropy is best modeled 
with a pair of Gaussians, each centerd ~ 50° away form 
the heliotail direction (see ^ in Table I). The amplitude 
of the southern excess (6i) is as large as ~0.2 % more 
than twice the amplitude of the northern excess (62), 
possibly indicating a significant north-south asymmetry 
of the heliosphere. The Tibet AS experiment succeeded 
for the first time to reveal a clear signature of the 
asymmetric heliosphere in the sidereal anisotropy of the 
multi-TeV GCR intensity. 
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